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ABSTRACT 

We describe an efficient algorithm for calculating the statistics of weak lensing by large-scale structure 
based on a tiled set of independent particle-mesh N-body simulations which telescope in resolution 
along the line-of-sight. This efficiency allows us to predict not only the mean properties of lensing 
observables such as the power spectrum, skewness and kurtosis of the convergence, but also their sampling 
errors for finite fields of view, which are themselves crucial for assessing the cosmological significance of 
observations. We find that the nongaussianity of the distribution substantially increases the sampling 
errors for the skewness and kurtosis in the several to tens of arcminutes regime, whereas those for the 
power spectrum are only fractionally increased even out to wavenumbers where shot noise from the 
intrinsic cllipticities of the galaxies will likely dominate the errors. 

Subject headings: cosmology: theory - gravitational lensing - large-scale structure of universe 



1. INTRODUCTION 

Weak lensing of background galaxies by foreground 
large-scale structure offers an opportunity to directly probe 
the mass distribution on large scales over a wide range of 
redshifts. In this paper we describe an N-body based algo- 
rithm optimized for weak lensing calculations which can be 
run on workstation-class computers. The method is fast 
and efficient, allowing the exploration of parameter space 
and production of many realizations of a given model to 
assess the statistical significance of any result. 

Weak lensing of distant galaxies by large scale struc- 
ture shears and magni fies th eir images. As first p ointe d 



out by Blandford et al. ( |l99lj[ ) and Miralda-Escude ( |l991h , 
these effects are of order a few percent in adiabatic cold 
dark matter models making their observation challenging 
but feasible. Early predictions for the power spectr um o f 
the shear and convergence were made by Kaiser (1992) 
on the basis of linear perturbation theory. Likewise the 
skewness of the convergence in perturbation theo ry wa s 
computed by Bemardeau, van Waerbeke & Mellier (1997). 
Jain & Seljak (1997) estimated the effect of non-linearities 
in the de nsity through analytic fitting formulae (Peacock 
& Dodds 1996) and showed they substantially increase the 
power in the convergence below the degree scale. 

On subdegree scales, a full description of weak lensing 
therefore requires numerical simulations, the most natu- 
ral being N-body simulations. N-body codes are ideally 
suited for weak lensing calculations since on the relevant 
scales only gravity is involved, bypassing the need for a 
treatment of hydrodynamic and radiative transfer effects. 
The evolution of density perturbations into the non-linear 
regime by N-body techniques is now a well-developed field. 
The particle-mesh (PM) N-body technique provides an ef- 
ficient means of simulating the evolution of structure. Its 
speed makes it ideal for the rapid exploration of cosmo- 
logical models and the calculation of statistical properties 
of the lensing observables, e.g. the sampling variance on 
estimators of the power spectrum, skewness and kurtosis 
of the convergence. While Lagrangian perturbation the- 
ory is arg uably even more efficient (Waerbeke, Bernardeau 
& Mellier 1998), without the proper non- linear dynamics 



one cannot guarantee that the statistics are faithfully re- 
produced. 

The main drawback of the PM technique is the lack of 
angular dynamic range, due partially to the broad kernel 
that describes the efficiency with which structures along 
the l ine-of-sight lens the sources (Jain, Seljak & White 
1999). We show here that this problem may be in large 



part overcome by tiling the line-of-sight with simulations 
of increasing resolution. 

The lensing signal is calculate d by r ay-tracing through 
the simulat ions 
& Ostriker Il99 



Blandford et al. 1991, Wambsganss, Cen 
Jain et al. 1999; Couch man, Barber & 



Thomas |1998| ; Fluke Web ster~gHVlortlock |1998| Hamana, 
Martel & Futamase 1999). In the weak lensing regime, a 



key simplification is that one can use unperturbed pho- 
ton paths to perform the relevant line-of-sight integrals, 
elimi nating the need for explicit ray tracing (Blandford 



et al. 1991; Hui, p rivate communication; see e.g. Bartel- 
mann & Schneider |2000 for a discussion) . This allows one 
to incorporate the lensing right into the time evolution of 
the code, eliminating the need to output the density field 
along the way and allowing very dense sampling of the in- 
tegrals. While the evaluation of the convergence along an 
unperturbed path is self-consistent within the framework 
of weak lensing, the results must be checked against a full 
ray traci ng method. The simulations reported in Jain et 
al. (1999) suggest that the approximation is good to 10 -3 
in power for models similar to the one reported in this 
paper. 

In this paper we shall concentrate on a specific cosmolog- 
ical model. It is a cosmological constant cold dark matter 
model (ACDM) with fl m ~ 0.3 = 1 — ft\, a scale-invariant 
spectrum of adiabatic perturbations (n = 1) with a mat- 
ter power spect rum described by the fitting function of 
Bardeen et al. ( |1986| ) with r BB KS = 0.2. The model is 
normalized to the COB E 4-year data using the method 
of Bunn & White ( |1997| ). This corresponds to cr 8 = 1.2, 
slightly above that infer red from the ab unda nce of rich 
clusters (Eke et al. |1998| , Viana & Liddle |1999| ). 

The outline of the paper is as follows. In we describe 
our implementation of a PM code and lensing evaluation. 
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In §[|we introduce the tiling technique. We present results 
for the power spectrum of the convergence and sampling 
errors in its estimation in §^ and analogous results for the 
skewness and kurtosis of the convergence §|[ A comparison 
of our tiling method and those based on single simulations 
is presented in We conclude in §[z| 

2. THE PM-LENSING CODE 

To evolve the dark matter distribution in the non-linear 
regime, we use a particl e-mes h (PM) code descr ibed in 
detail in Meiksin et al. ( [L999D and White ( [19991) . The 
simulations reported here use either 128 3 or 256 3 particles 
and a 256 3 or 512 3 force "mesh". The initial conditions 
are generated by displacing particles from a regular grid 
using the Zel'dovich approximation. The simulations are 
started at 1 + z = 35 and evolved to the present (z = 0) 
using adaptive steps in the log of the scale factor a = 
(1 + z)~ x . The force on each particle is calculated from 
the density using Fourier Transform (FT) techniques with 

a kernel —k/k 2 . The gridded fields are computed from the 
particle d ata u sing CIC charge assignment (Hockney & 
Eastwood 1981). The time step is dynamically chosen as 
a small fraction of the inverse square root of the maximum 
acceleration, with an upper limit of Aa/a = 3 per cent per 
step. The code typically takes 200-300 time steps between 
1 + z = 35 and z = 0. 

The new ingredient, beyond simple N-body evolution 
of the density field, is the calculation of the convergence 
along a bundle of rays. In the weak lensing approximation, 
calculation of this scalar quantity in any field is sufficient 
to enable calculation of all of the other quantities (e.g. the 
shear 7). We assume here for simplicity that the sources 
all lie at one redshift z s = 1. The code as written allows 
multiple source redshifts, but we restrict ourselves to the 
single source plane in this paper. 

Before the N-body evolution begins, we generate geo- 
desies, in code coordinates, for iVi os = 128 2 or 256 2 lines- 
of-sight by integrating 



dD 



1 



da 



a 2 H(a) 



(1) 



where D11 is the comoving distance parallel to the line-of- 
sight. The lines-of-sight originate in a square lattice at z s 
and converge upon an observer situated at the center of 
one face of the box at z = 0. We further demand that the 
field of view never subtend more than a box length to avoid 
introducing artifacts due to periodic boundary conditions. 
We make the small angle approximation and assume that 
-D|| lies parallel to the z-axis for all rays, thus the coordi- 
nates perpendicular to the line-of-sight scale linearly^] with 

The N-body code is then run, and once the evolution 
reaches a redshift of z Sl we integrate the lensing equation 



k(x ± ) = D s I dD\\ f(l - t)V\<$>(x) 



(2) 



in addition to the gravitational force. Here t = D{a)/D s 6 
[0, 1] is the dimensionless distance to the source. For mul- 

lr This is appropriate for the flat universes we deal with in this 
paper. In a curved universe, the angular diameter distance must be 
used to calculate the "opening distance" of any ray from the center 
of the box as a function of redshift. 



tiple sources one can replace the kernel for source i with 
t(l — t) with t(t S i — t)/t S i where t s = D S i/D s is the distance 
to the zth source in units of the distance to the furthest 
source D s . 

The source V 2 L < i ) (x) is calculated in the box using FT 
methods under the small-angle approximation. The par- 
ticles are assigned to the n earest point on a grid (NGP; 
Hockney & Eastwood 198l| ) to obtain the density distri- 
bution. The FT of this distribution, 5k, is then multi- 



plied by -|fi r 



|fi m i/Qa _1 A: 2 /k 2 and the transform inverted. 
Within each time step, we assume that the potential is 
slowly varying <i>(a 4- 5a) ss &(a). Since time steps are 
separated by Aa/a ~ 0.01, much less than the expansion 
time on which $ varies, this is a very good approxima- 
tion. The integral is evaluated by taking N points along 
each line-of-sight and time step assuming the potential is 
frozen. By increasing N, we find that N ~ 10 2 dynami- 
cally chosen points suffices for convergence. This sub-step 
integration range runs from the a of the last time step in 
the code to the present a. The integral in Eq. (||) is there- 
fore densely sampled and the n correctly evaluated at the 
a corresponding to the box redshift. 

We first test the Limber approximation (see Appendix) 
which says that only modes perpendicular to the line-of- 
sight contribute to the integral in Eq. (J2j> . In this approxi- 
mation, the 2 dimensional Laplacian can be replaced with 
a 3 dimensional Laplacian which in turn can be expressed 
in terms of the density perturbation through the Poisson 
equation: 



k(x x ) = ^n m H 2 D s J dD\\ t(l 



t)[5{x)/a]. (3) 



Using integration by parts the error induced by th is re- 
placement should be £>($) ~ 10 -5 (Jain, et al. |l999| ). 

We have run a 256 3 PM simulation of our ACDM model 
in a 125/i _1 Mpc box using Eqs. (||) and (||) to compute k 
in a grid of 256 2 lines-of-sight. The two track each other 
very well. The power spectra computed from the two fields 
are almost identical, as are the histograms of k. In a line- 
of-sight by line-of-sight comparison Eqs. (|J) and (||) return 



Table 1 



Tiling Solution 



Gout 


Lbox 


flout 


Lbox 


0.537 


245 


0.822 


75 


0.577 


245 


0.841 


75 


0.610 


195 


0.861 


75 


0.646 


195 


0.881 


75 


0.675 


155 


0.902 


75 


0.707 


155 


0.924 


75 


0.732 


120 


0.946 


75 


0.759 


120 


0.970 


75 


0.780 


95 


0.994 


75 


0.803 


95 


1.000 


75 



NOTES. — The tiling solution for our ACDM model assuming a single 
source redshift z B = 1, i.e. a B = 0.5, that uses 6 box sizes and 20 tiles. 
The column a ou t gives the scale-factor at which each tile is output. 
(A tile contains that part of the integration of re lying between the 
last output and a ou t.) The size of the simulation box used for that 
tile (in h~ 1 Mpc) is also given. 
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values for k that deviate by at most 0.03 and on average 
(rms) 0.003. For comparison, the rms fluctuation on the 
grid scale in these planes is nearly an order of magnitude 
larger than this: a K ~ 0.02. Some of this scatter is no 
doubt induced by our small-angle approximation in com- 
puting Vj_, while some comes from the finite size of the 
box. Since the integration has traced across the box 19 
times and at each edge we can pick up a term C($), this 
level of variance agrees roughly with our expectations. In 
the absence of finite box size effects, we expect Eqs. (|2j) 
and (^) would match even more closely and there is rea- 
son to believe that our eva luatio n of Eq. (||) is the more 
accurate (see also Stebbins 1999). 

Since Eq. (||) is less computationally expensive, we will 
adopt it from this point on. The fact that this approxima- 
tion works well shows that the integral Eq. (||) is sensitive 
only to those modes in the box which are perpendicular 
to the line-of-sight. This is an important point to remem- 
ber when considering questions of sample- or run-to-run 
variance. 

Finally the entire bundle of rays is rotated at a random 
angle to the box faces and placed at a random offset from 
the box center. This ensures that the rays do not trace par- 
allel to the edges of the simulation box and the grid used 
to define the density. Since the simulation uses periodic 
boundary conditions, we actually compute the density in 
a periodic universe. Thus each time a ray leaves the box 
it is remapped into it using periodicity. 

3. TILING 

A photon from z ~ 1 traverses many Gpc on its way to 
us whereas the large-scale structure responsible for lensing 
spans the Mpc range and below. Simulating the full range 
of scales implied is currently a practical impossibility. A 
solution commonly employed in the literature is to recycle 
the output of a single smaller simulation, i.e. sum the 
contributions of the density, projected to the midplane, of 
the given simulation at a series of discrete redshifts. We 
propose here a "tiling" alternative that addresses three 
potential problems with the traditional technique: the lack 
of statistical independence of the fluctuations, the loss of 
angular resolution in the projection, and the discreteness 
of the projection. 

We maintain the the statistical independence of the fluc- 
tuations by employing multiple independent simulations to 
tile the line-of-sight. We are free then to adjust the sizes 



Table 2 



Number of Simulations 



^box 


A^256 


N b l2 


245 


76 


10 


195 


20 


10 


155 


20 


10 


120 


20 


10 


95 


21 


10 


75 


36 


15 



NOTES. — The sizes of the simulation boxes used (in h~ 1 Mpc) and 
the number of independent boxes of resolution of that size, with both 
256 3 and 512 3 mesh resolutions. 



of the simulation boxes and in particular can make them 
smaller and smaller as the rays converge on the observer 
(see Table ||). Recall that as the lines-of-sight converge, 
they probe ever smaller physical separations for a given 
angular separation. Since the lensing kernel is so broad, 
even structure quite close to the observer contributes to 
the signal. In fact, the non-linear amplification of struc- 
ture at low-z implies that on large angular scales the lens- 
ing kernel is skewed toward the observer (see Fig. 1). 

Specifically for each simulation box, the code outputs 
the contribution to the k planes at specified redshifts (see 
Table |l|), typically spaced in a so that the photons tra- 
verse the box once between each output. Note that this is 
not the same as simply computing the projected density 
at the midplane. The full integral, with the evolution of 
the potential and the geometry of the rays etc, is being 
computed within each tile. After each output the offset 
and random orientation of the rays are chosen anew and 
the integration is started afresh for the next segment. 

If we simulate only a single box, the integral of Eq. (j|) 
is simply given by the sum of the planes from that box. 
However with multiple simulation boxes run with the same 
tiling scheme, we can then construct our final k plane as 
the sum of planes from different simulations. In practice, 
we shrink the box size so that it fits in the field-of-view at 
the cndplanc until it reaches the non-linear scale. The non- 
linear scale must be within the box at the relevant epoch 
to ensure that the PM code evolves the density correctly. 
Nonetheless we shall demonstrate that this box resizing 
technique is very effective by comparing results from a 
series of low resolution (256 3 ) simulations to our higher 
resolution (512 3 ) simulations done in a box of a single size 
(§§■ 

The result at the end of the simulation(s) is a grid k 
along lines-of-sight spaced equally in angle. We then cal- 
culate the shear from this grid by using 
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Fig. 1. — The contribution _tp as a function of scale-factor for 
i = 100, 300, 1000 from Eq. (§). 



Fig. 2. — (left) An image, 6° on a side, of the convergence ft from a single realization of our tiling solution. The grey scale is linear from 
ft = —0.05 to 0.15. (right) The shear field, 74, derived from the left panel. The lines have been exaggerated, the amplitude of the shear is at 
the percent level as in (a). 



72 



2lil 2 



(5) 



where k is the 2D FT of the convergence field,^ and £ = 
(£1,(2) is the Fourier variable conjugate to the position on 
the sky. 

We show in Fig. || the convergence n and the derived 
shear field 7j, from one of the 512 3 simulations using the 
tiling scheme described in Table [|. The field is 6° on a side 
and contains 256 2 lines-of-sight. From our multiple simu- 
lations, we are able to generate many independent fields 
of this size and resolution. In the following sections we 
discuss the statistics of these fields based on 512 random 
combinations of the tiles listed in Table || for both the low 
and high resolutions sets. 

4. POWER SPECTRUM 

Fig. |^ shows the angular power spectrum of ft, com- 
puted from the tiling simulations, as compared to the semi- 
analytic result (see Appendix), 



9tt 

4£ 



si m H*D*] 2 J^t 3 (i-tf 
(k = i/D {l , a y 



A 



(6) 



where A^ ass (fc) = fc 3 P(fc)/(27r 2 ) is the contribution to the 
mass variance per logarithmic interval physical wavenum- 
ber and analogously A 2 (£) = 1 2 Ci/{2-k) is the contribution 



to fCj ms per logarithmic interval in angular wavenumber 

2 Because the field is not periodic, it is important to zero-pad the 
FT array before computing ft. 




100 



1000 



Fig. 3. — (top) The angular power spectrum, £ 2 Ce/(2ir) or A 2 , 
vs. multipole number I, for the convergence ft from our tiling sim- 
ulations. We also show the semi-analytic prediction from Eq. (H) 
using both linear theory and the non-linear power spectrum. Tne 
shot-noise contribution assuming n = 2 X 10 5 galaxies per square de- 
gree each with an rms ellipticity 7 rm s = 0.4 is also shown, (bottom) 
The patio of the simulation results to the (non-linear) prediction of 
Eq. (H). 
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Table 3 



Simulation Properties 



^mid 


Weight 


Lbox 


#box 


-^mesh 


^mcsh 


^part 






(/i _1 Mpc) 


(arcmin) 


/ 7 1 i \ 

(h i kpc) 


(arcmin) 


(1O 9 M ) 


0.518 


0.05 


245 


385 


479 


0.75 


73 


0.557 


0.13 


245 


433 


479 


0.85 


73 


0.593 


0.19 


195 


389 


381 


0.76 


37 


0.628 


0.22 


195 


438 


381 


0.86 


37 


0.660 


0.24 


155 


393 


303 


0.77 


19 


0.691 


0.25 


155 


444 


303 


0.87 


19 


0.719 


0.25 


120 


388 


234 


0.76 


8.6 


0.745 


0.25 


120 


438 


234 


0.85 


8.6 


0.769 


0.23 


95 


391 


186 


0.76 


4.3 


0.792 


0.22 


95 


441 


186 


0.86 


4.3 


0.81z 


O.zU 


(a 


on a 
394 


A AC 

146 


U.77 


Z.l 


0.831 


0.19 


75 


444 


146 


0.87 


2.1 


0.851 


0.17 


75 


510 


146 


1.00 


2.1 


0.871 


0.15 


75 


599 


146 


1.17 


2.1 


0.891 


0.13 


75 


726 


146 


1.42 


2.1 


0.913 


0.11 


75 


920 


146 


1.80 


2.1 


0.935 


0.08 


75 


1257 


146 


2.45 


2.1 


0.958 


0.05 


75 


1981 


146 


3.87 


2.1 


0.982 


0.02 


75 


4673 


146 


9.13 


2.1 


0.988 


0.02 


75 


6875 


146 


13.43 


2.1 



NOTES. — As a function of the scale factor at the middle of each tile: the weight, i(l — t), at the tile midpoint, the box size used for that tile 
and the size of the force mesh, the angular size of the box and mesh and the particle mass. These numbers are for the 512 3 simulations. For 
the 256 3 simulations the mesh size should be doubled and the mass per particle increased by 2 3 . 



(or equivalently multipole) t. We also show, in Fig. |, the 
power spectrum from the first 5 realizations, to emphasize 
the scatter from field-to-field. 

In evaluating Eq. (Q) , we use the method of Peacock & 



< 



io- 




100 



1000 



Fig. 4. — As Fig. |3J, but with the power spectrum from 5 different 
realizations shown to emphasize the scatter_Trom field-to-field. The 
thick solid line is the prediction from Eq. <M). 




0.1 1 
k (h/Mpc) 



Fig. 5. — (top) The 3D mass power spectrum from our ensemble of 
simulations as compared to the fitting function of Peacock & Dodds 
(1996). Here and only here the error bars represents the error on 
the mean computed from averaging over our many realizations of 
the model, (bottom) The ratio of the N-body results to the fitting 
formula. 
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Dodds (1996) to compute the non-linear power spectrum 
as a lunction of scale-factor. Comparison with the average 
power spectrum from our simulations (e.g. Fig. ^| at z = 0) 
shows agreement at the 10% level for the range of redshifts 
and scales resolved by simulation (k <; 5/iMpc _1 ). 

The loss of power on large scales (small £) is a result 
of our FT based analysis routines and the 6° x 6° field 
of view. To test this we generated gaussian fields with 
the angular power spectrum of Eq. (^) and with much 
larger areas. When analyzing 6° x 6° subfields the same 
\ow-£ suppression as in Fig. |3J was seen and comes from 
"windowing" the map by the field of view. 

The roll-off at high-i? in Fig. |^ is as expected from the 
resolution of the N-body code. The PM code resolves 
scales (in k) down to approximately &Nyquist/3 with a slight 
dependence on the spectral index of the model. The small- 
est k we can simulate is 27r/Lbox so in a 512 3 simulation 
we would expect a dynamic range in k of 256/3 ~ 90. 
The projection from physical scale to angular scale is not 
unique but rather has a finite width "kernel" (see Eq. |J). 
In our case the width is roughly a factor of 3 in scale. So 
to fully resolve a given £-mode, we need to resolve a factor 
of 3 higher in physical scale than that mode projects at 
the mid-plane. This further reduces our dynamic range to 
a factor of 30. Because we cannot arbitrarily reduce the 
box size without the fundamental mode going non-linear 
by the present, our tiling is inefficient at low-z and our 
actual dynamic range is closer to a factor of 20-25, as can 
be seen in Fig. |[ 

The error bars in Fig. |^ are the sampling errors for an in- 
dividual 6° x 6° field of view as estimated from the scatter 
of the full suite of simulations. This should not be confused 
with the much smaller error on the mean power spectrum 
of the suite. Sampling errors for different survey dimen- 
sions scale roughly as the ratio of the dimensions and the 



variance as the ratio of the survey areas. 

As demonstrated in Figs. [| g, even though sampling- 
errors only fully converge to that of a gaussian random 
field with the same power spectrum for I <; 300, the non- 
gaussian contribution to the errors remains in the few tens 
of percent out until at least I <; 1000 (in qualitative ac- 
cord wi th an alytic estimates, see Scoccimarro, Zaldarriaga 
& Hui 1999). We have checked that the deviations from 
gaussianity are only weakly dependent on the binning cho- 
sen for this range in I. This can also be seen be examining 
the covariance of the binned power spectrum estimators 
shown in Table || As with the variance, the covariance 
deviates from the gaussian limit beginning at I ~ 300 and 
grows at a moderate rate through I ~ 1000. The bins are 
correlated even in the gaussian limit by the limited field of 
view: the fundamental mode implies a spacing of A£ — 60. 

The full distribution of the power spectrum estimator 
also becomes moderately less well characterized by its vari- 
ance for I > 300. In Fig. |^, we show the histogram of val- 
ues from the simulations. The probability of outliers on 
the low side decreases due to the skewness of the distri- 
bution whereas that on the high side remains reasonably 
well characterized by the variance for 2 and 3c outliers. 
These probabilities are with respect to a Gaussian sky ar- 
tificially set to the same variance for the power spectrum 
estimator. This should not be confused with the expec- 
tations from the Gaussian prediction for the variance: for 
the I — 738 bin, a > 2a deviation from the mean power 
with respect to the Gaussian standard deviation occurs in 
one quarter of our tiles. 

Beyond I ~ 1000, the nongaussian contributions to the 
variance, covariance, and tails of the distribution of the 
power spectrum estimators becomes substantial. However, 
this is also the point at which shot noise from the intrin- 
sic ellipticity of the galaxies begins to exceed the sample 
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Fig. 6. — (top) The standard deviation SA'f. in our binned es- 
timates of as a function of I. The solid line is from our 512 3 
simulations, dotted is our 256 3 simulations and the short-dashed 
line is from Gaussian fake skies with_the same power spectrum. The 
choice of binning is given in Table W. (bottom) The ratio of errors 
with shot-noise to pure sampling errors for the three cases above. 



Fig. 7. — The histogram of Af.(£) for two different bins from our 
higher (solid) and lower (dashed) resolution simulations, and fake 
fields generated with Gaussian statistics (dotted). 
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variance. The shot noise power spectrum is (Kaiser 1998) 

7 2 

[ rms /rj\ 

noise _ : \ ' J 



where n is the number density of the sources and 7 rms is 
the rms intrinsic shear in each component. The shot noise 
spectrum for n — 2 x 10 5 deg -2 and 7 rms = 0.4 is shown in 
Fig. 0. The noise bias in the measurements of the power 
spectrum can be subtracted off at the expense of increasing 
the variance of the estimator for each l-mo&e 



SC 



| total 



= SC, 



e Ik 



4C,C n 



2Ci 



(8) 



For our binned estimators, the sample variance is reduced 
by \f~N statistics so that the total fractional variance is 



SCt 
Ct 



total 



5Ct 



1 



2n2 



NtC 



(iC e C noisc + 2C 2 oise ) 



(9) 

where the first term is the result from our simulations 
(without shot-noise) and the sum in the second term is 
over the N/> independent modes in the bin. The num- 
ber of independent modes for a given £ is approximately 
(2£ + l)/sk y , where / s k y is the fraction of sky covered by 
the field of view (/ s ky ~ 10 -3 for our fields). We show 
the effect of shot noise on the sample variance in Fig. |[ 
We have tested these approximations against monte carlo 
realizations of the shot noise and found good agreement. 

The combination of these results imply that techniques 
based on gaussian assumptions for power spectrum esti- 
mation are fair approximations at lea st in t he con text of 
this AC DM m odel (e.g. Kaiser |1998| Seljak |1998| ; Hu & 
Tegmark |l999| ). 



5. SKEWNESS AND KURTOSIS 

Figure ^ shows the co-added histogram of n, smoothed 
on 5' and 10', from 512 tiling solutions. The non-gaussian 
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Fig. 8. — The histogram of re, smoothed with a top-hat filter of 
radius 5' and 10'. The solid lines are from our 512 3 simulations 
while the dashed lines are from our 256 3 simulations. 



nature of the distribution is apparent in this figure, as is 
the 1ow-k cutoff enforced by <5 mass > — 1. The large number 
of tiling solutions we have run allows us to probe the dis- 
tribution well into the tails. Clearly the higher and lower 
resolution simulations agree well on these scales. Our abil- 
ity to simulate many k planes allows us to study the statis- 
tics of the moments of this distribution. In this section, 
we examine the lowest order moments beyond the 2-point 
function: the skewness and kurtosis. 

5.1. Simulation Results 

From the two dimensional angular grid of the conver- 
gence k, we calculate the skewness and kurtosis on an an- 
gular scale cr. We first smooth the grid with a pixelized 
tophat window W a with FT techniques 



(10) 



and eliminate edge effects by zero padding the array and 
discarding the data that is convolved with the zero padded 
region. We then calculate the skewness 



\2 ' 



and the kurtosis 



(ii) 



(12) 



for two different averaging schemes: averaging over pixels 
in a given 6° x 6° field and averaging the pixels over all 
fields. 

As can be seen in Figs. || and [H)[ even a 6° x 6° field 
suffers from large sample variance on scales of tens of ar- 
cminutes. Like the power spectrum estimators, we expect 
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Fig. 9. — The skewness, 53, as a function of (top hat) smoothing 
scale. The squares are results from the 512 3 simulations while the 
circles are from the 256 3 simulations. Filled symbols indicate the 
skewness computed from the set of generated re planes while open 
symbols with error bars indicate the mean and variance of S3 for 
each plane. The points are offset slightly for clarity. The solid line is 
a semi-analytic estimate (L. Hui, private communication), discussed 
in more detail in the text. 
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Table 4 



Power Spectrum Covariance 
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138 


(0.23) 


1.00 


0.31 


0.21 


0.09 


0.14 


0.16 


0.18 


0.15 


0.22 


194 


(0.04) 


(0.22) 


1.00 


0.26 


0.24 


0.28 


0.17 


0.15 


0.19 


0.16 


271 


(-0.02) 


(-0.03) 


(0.17) 


1.00 


0.38 


0.33 


0.34 


0.27 


0.19 


0.32 


378 


(-0.01) 


(0.02) 


(0.04) 


(o.n) 


1.00 


0.45 


0.38 


0.33 


0.32 


0.27 


529 


(0.01) 


(-0.01) 


(-0.07) 


(-0.02) 


(0.02) 


1.00 


0.50 


0.48 


0.36 


0.46 


739 


(0.04) 


(-0.03) 


(-0.01) 


(-0.02) 


(-0.02) 


(0.13) 


1.00 


0.54 


0.53 


0.50 


1031 


(0.07) 


(0.01) 


(0.07) 


(-0.03) 


(0.03) 


(0.08) 


(0.04) 


1.00 


0.57 


0.61 


1440 


(-0.03) 


(0.02) 


(0.04) 


(-0.04) 


(0.05) 


(-0.07) 


(-0.04) 


(-0.03) 


1.00 


0.65 


2012 


(-0.02) 


(-0.04) 


(0.03) 


(0.03) 


(0.03) 


(-0.02) 


(0.03) 


(-0.07) 


(0.02) 


1.00 



NOTES. — Covariance of the binned power spectrum estimators. Upper triangle displays the covariance found in the 512 tilings of the 512 3 
simulations. Lower triangle (parenthetical numbers) displays the covariance found in an equal number of gaussian realizations. The finite 
6° X 6° field of view couples the power spectrum estimators over At ~ 60 in both cases, whereas non-linear dynamics couples the estimators 
in the simulations at high I. 



the sample variance to scale roughly with the survey area. 
Through generating Gaussian fields with the same power 
spectrum, we find that the sampling errors for S3 and S4 
are a factor of 2 and 7 larger than the Gaussian limit re- 
spectively for a ~ 10'. 

The difference in the moments computed from averaging 
Sn in each field compared to computing Sn usin g all o f 
the field has been stressed by Hui & Gaztanaga ( |1999| ). 
The bias increases as the field-to-field variance increases 
as can be seen by comparing large and small smoothing 
scales in Figs. and |l0|. (In our simulations we found 
that the value 01 Sn computed using the moments of all 
the fields fluctuated more with increasing numbers of runs 
than the mean of the Sn computed from moments within 
each field.) This large sampling errors should be borne 
in mind when employing S3 measurements to distinguish 
between cosmological model. 

Comparison of the 512 3 simulations with the 256 3 sim- 
ulations indicates that the N-body calculation has con- 
verged on a scale of 10', both in the moments themselves 
and in the sampling errors. The two sets of simulations 
begin to diverge in their fractional standard deviation near 
5', suggesting that the higher resolution simulations may 
even be reliable down to 2.5'. Fig. |Il] shows the diver- 
gence between the simulations in S3 is m the high S3 tail, 
which may be due to resolution or may indicate too few 
higher resolution simulations have been run. We have also 
checked that the 75/i _1 Mpc are large enough to provide an 
adequate sample of the non-linear scale for these purposes. 
Omiting these simulations and completing the tiling with 
95/i~ 1 Mpc simulations produces a negligible change in S3 
at 10'. 

As Fig. O shows, the kurtosis increases above the 
(re 2 )S4 = 3 Delow 10'. As this is the number expected 
for (k*) / (k 2 ) for a gaussian field, it marks the regime 
where the distribution becomes significantly non-gaussian 
in the 4th moments. However we detect no similar dra- 
matic rise in the power spectrum errors at I ~ 1000 (§0). 

Finally we have simulated the effect of shot-noise on the 
variance of S3 and S4. In the presence of shot-noise we 
define estimators of Sat in analogy with Eqs. (O, O) but 



which subtract the contribution of the shot-noise to (k™). 
For example if k! is the measured value of K a including 
shot-noise with variance (e„), we define 



./ 3 



« 2 > 



(13) 



Using these estimators and adding simulated shot-noise 
to our planes we find that the estimators are unbiased 
and their standard deviations are only slightly increased 
(< 16% for S3 and <; 6% for S4) even on scales as small as 
2.5'. This is not too surprising since with 2 x 10 5 galaxies 
per square degree the shot-noise power only surpasses the 
signal power in our model on scales smaller than 1.3', and 
we have shown that the sample variance on S3 and S4 is 
significantly enhanced by the non-Gaussianity of the dis- 
tribution. Artificially increasing the noise by a factor of 4 
does lead to an increase in the variance of £3 and 54, but 
the estimators remain unbiased. 

5.2. Comparison to Previous Results 

These results make sense physically, but it would be 
useful to compare with previous work. On the scales we 
are working perturbation theory is not adequate, so the 
best comparison is with o ther s imulations, the closest be- 
ing those of Jain, et al. (199£). Unfortunately a direct 
comparison with their work is difficult. While our model 
differs slightly from theirs, we have run a smaller set of 
simulations of their exact model and find that S3 are not 
strongly affected by the slight changes. However we do not 
have the dynamic range to reliably estimate the skewness 
on 1' scales, and their 3.5° x 3.5° field is sufficiently small 
that they have large sample variance on 10' scales. Using 
our analysis software on one field from their simulations 
(B. Jain, private communication), our skewness is approx- 
imately 20% lower at 5' than theirs. We compare at 5', 
which is the edge of our reliable range, because the sample 
variance from their small fields makes comparison difficult 
above this scale. Indeed, in the plane we have analyzed 
their skewness peaks at 10' before dropping precipitously. 
In the distribution of S3 in our higher resolution (512 3 
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mesh) simulations (see Fig. [Tl]) only 5% of our plan es hav e 
5*3 as high or higher than the plane from Jain et al. ( 1999| ). 
Crudely accounting for the increased sample variance due 
to their smaller field by scaling the distribution by 6°/3.5°, 
raises this number to 10 — 15%. While this is not highly 
improbable, the difference may still be du e to systematic 
differences in the codes. Jain, et al. ( 1999 ) also performed 
some PM runs in a 64/i _1 Mpc box with a 256 3 force mesh 
and found results ~ 20% lower than their P 3 M results at 
5' (Seljak, private communication). Whether this discrep- 
ancy is due to the small box size they used, syste matic 



difference between P M and P 3 M (e.g. Splinter et al. [1998 
Jain & Bertschingei n|l99 <j ) or sample variance is not clear 
We show in Figs. |9|, |11| the prediction of a semi-analytic 
calculation (L. Hui, private communication) based on 
hyper-exte nded perturbation theory (HEPT; Scoccimarro 
& Frieman 1999). The agreement is at the level one would 
expect from the approximation used. To check this we 
have calculated the skewness and kurtosis of the density 
field (at z — 0.4, the peak of the curves in Fig. 1) in 
our 155/i _1 Mpc boxes with 256 3 particles and 512 3 force 
mesh. For each of the 10 simulations, we binned the par- 
ticles onto a 512 3 grid using NGP assignment (the results 
do not depend on this choice), then smoothed this grid 
using a 3D analogue of Eq. (10) with a top hat radius R. 
The moments were computed by averaging powers S over 
the 512 3 grid sites. Again we computed the average Sn 
over the simulations and the "global" Sn from combining 
the moments from all the simulations: 



S 3 (R) 



S±(R) 



(14) 
(15) 



3 The discrepancy noted by these authors for the n = —2 spectrum 
is not directly relevant here, since n is less negative on the scales of 



interest. We have shown this specifically in Fig. |5| The general point 
remains valid however. 



Our results are shown in Fig. |l2] as a function of radius, 
along with the variance in the density. Also plotted (dot- 
ted) are the predictions of HEP T as used in the semi- 
analytic lensing calculation (Hui 1999| ) and t he variance 
(dashed) predicted by Peacock & Dodds ( |1996| ). The non- 
gaussianity in the lensing signal may be generated at lower 
z than the peak in Fig. 1, so we have also calculated S3 and 
S4 from our z = data. The results are consistent with 
little or no evolution in S3 and S4 since z = 0.4, tho ugh th e 
variance grows as predicted by Peacock & Dodds (1996). 

The level of disagreement is sufficient to explain the dis- 
crepancy in Fig. 0. The stated realm of validity for the 
HEPT result is for variances > 100 and indeed our results 
for S3 and S4 of the density in this regime agree better 
with the prediction (see Fig. 12). For lower variances, one 
expects both moments to be smaller and the power law 
approximation to the mass power spectrum inherent in 
HEPT to break down. Preliminary results from a treat- 
ment that includes these two effects (R. Scoccimarro, pri- 
vate communication) indeed agree better with our lensing 
results: S3 « 115 at 5', compared with our 111, with a 
gradual decline to S3 ~ 80 at 70' — 80' before an increase 
bac k to t he perturbation theory results of Bernardeau et 
al. (1997). In any case, the difference shown in Fig. |lj 
is very likely the cause of the discrepancy with the semi- 
analytic calculation. 

6. TILING TESTS 

Due to our "tiling" method, and the large number of 
simulations that we have run (Table |), we are able to 
systematically examine the dependence of our various re- 
sults on the volume of space sampled by the simulations. 
Of particular interest is the following question: how is 
the sample variance associated with each field affected by 
tracing repeatedly through a single simulation? We have 
attempted to answer that question for various statistics 
with our large ensemble of simulations. 




Fig. 10. — As Fig. y but for the kurtosis, 54. The kurtosis has 
been scaled by (k 2 ) for display purposes. 
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Fig. 11. — The histogram of Ss(5') values from our higher (solid) 
and lower (dashed) resolution simulations. Also shown are the pre- 
dictions from_HEPT (see text) and one plane of the simulations of 
Jain et al. (1999). 
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We first looked at the statistics of the power spectrum. 
Using our 76 large boxes (245ft.~ 1 Mpc at 256 3 force reso- 
lution) we checked that the mean, variance, and pdf of the 
power spectrum (for 4 different binnings) was the same 
whether we shuffled the tiles between boxes or used each 
box in isolation. These large and lower-resolution boxes 
are not fully resolving the structure at late times, but this 
is not of great concern as the small scale structure that 
is missing is unlikely to be correlated over large scales as 
required to cause an effect in this test. 

For the power spectrum, repeatedly tracing through 
a 245/i _1 Mpc box provides the same distribution as our 
tiling method (though multiple boxes are still needed to 
assess sample variance). The same test for our smaller 
75/i~ 1 Mpc also shows no statistically significant effect. 
This is significant since as the box shrinks the volume of 
space sampled is reduced and sampling becomes a larger- 
issue . Likewise for S3 and S4 no significant difference be- 
tween repeated tracing and tiling was found. 

These results shed light on another possible concern: 
that the rays in these simulations trace through boxes 
which are joined "discontinuously" at their edges. Fig. || 
shows that this does not affect the mean value of the 2- 
point function reproduced by the code. Our multiple sim- 
ulations allow us to go further however. A comparison 
of the tiling simulations with the repeated tracings of the 
245/i~ 1 Mpc box allows us to test the effect of a different 
number of "matchings" along the line-of-sight. On the 
scales where both resolve the structure, we find conver- 
gence in the mean value of S3 (the values of S4 are too 
noisy to allow a strong statement) and the pdf of Aj; with 
4 different binnings. This is suggestive that box matching 
is not a major source of error, though we cannot test that 
this is true on very small scales due to lack of resolution 



in our simulations. 

Tiling does increase the angular resolution of our simu- 
lations. In Fig. |l3| we compare the angular power spectrum 
from our higher resolution (512 3 mesh) PM runs in boxes 
of size 245ft.~ 1 Mpc with our tiling solution at lower res- 
olution (256 3 mesh), using the shrinking box. One can 
see that allowing the box to shrink along the line-of-sight 
produces considerable gains in angular resolution. Unfor- 
tunately the need to keep the fundamental mode of the 
box linear at all times restricts the size of the low-z boxes 
and limits the gain in angular resolution which can be 
achieved by this method (to a factor between 2 and 3). 
While larger fields of view are easily simulated, the mini- 
mum size of the low-z boxes restricts the smallest angular 
scale that can be probed. This gain is enough, however, 
to make PM codes viable for a rapid exploration of pa- 
rameter space on workstation class machines (cf. Jain, et 
al. |1999|). 



7. DISCUSSION 

We have described an efficient algorithm for calculat- 
ing the statistics of weak lensing by large-scale structure 
in N-body simulations. By working with the unperturbed 
paths, our method is extremely simple to implement and 
can be done at the same time as the N-body run(s). This 
gives one the ability to simulate a large volume and sam- 
ple the line-of-sight integration densely in both space and 
time. Contrast this with more traditional ray tracing tech- 
niques which use only tens of lens planes and project the 
entire density distribution in the box onto a single lens 
plane for each time step. Neglecting the deflections is cer- 
tainly self-consistent within the weak lensing approxima- 
tion. Analytic arguments (Bernardeau et al. 1997) and 
explicit ray-tracing simulations (Jain, et al. 1999) furthcr- 
more imply that corrections due to deflections are small 
for our purposes. 
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Fig. 12. — The variance, skewness, S3, and kurtosis, 54, of the 
density field at z = 0.4 as a function of (top hat) smoothing scale. 
These results are from our 512 3 simulations, as described in the text. 
Filled symbols indicate the 5jv computed from 10 boxes, while open 
symbols with error bars indicate the mean and variance of 5jv for 
each box. The dotted lines arpjthe predictions based on HEPT, used 
to generate the eurye in-Fig 0. The dashed line is the prediction of 
Peacock & Dodds (1999). 
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Fig. 13. — Comparison of resolution from tiling compared with 
single size boxes. The open circles represent the power spectrum of 
fields produced by "tiling" with the 256 3 PM simulations. The filled 
squares are the power spectrum deduced from our higher resolution, 
512 3 simulation, but using only the largest box: 245/i — 1 Mpc. Notice 
that tiling wins back the extra factor of 2 in resolution. 
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As with other simulation based results, numerical reso- 
lution and dynamic range are a serious issue. In particular 
the effects of finite force resolution can be seen in our re- 
sults below 2.5'. For weak lensing there are also problems 
introduced by the periodicity of the simulation box, which 
limits the size of the field of view that one can probe in a 
given simulation, in our case to 6° x 6°. 

We have described a technique, which we call "tiling", 
which allows us to use results from multiple realizations 
of a given model, and to match the size of the simulation 
box to the converging ray bundle to increase angular reso- 
lution for a fixed physical resolution. By varying the tiling 
scheme, we also tested the effects of discontinuities from 
joining the boxes and repeatingly tracing through the same 
simulation. We found no significant effect from either. 

With our suite of simulations, we are able to predict not 
just the mean properties of the models, but also their sam- 
pling errors. This is extremely important in assessing the 
statistical significance of future measurements. We have 
shown that the non-gaussian contribution to the errors on 



the power spectrum remains small out to I 



though the distribution of convergence, k, is clearly non- 
gaussian at 10'. We have quantified the (large) sample 
variance in estimates of higher moments of the k distri- 
bution which have been suggested as tests of the energy 
contents of the universe. Even with 6° x 6° fields the errors 
on the moments arc totally dominated by sample variance 
on scales above a few arcminutes with galaxy densities 
achievable in current observations. 

Since sampling errors scale inversely with the dimen- 
sions of the survey, a field of view in the tens to hundreds 
of square degrees will be crucial for extracting cosmologi- 
cal information on large scale structure from weak lensing 
surveys, especially for the non-gaussian signatures of mod- 
els. Nevertheless, due to the growing number of instru- 
ments with wide fields o f view , for example MEGACAM 
at CFHT (Boulade et al. |1998| ) and the VLT Survey Tele- 
scope at th e European Southern Observatory (Arnaboldi 
et al. 1998| ), the prospects for weak lensing in the era of 
precision cosmology remain bright. 
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APPENDIX 

LIMBER'S EQUATION 

In this appendix we provide a simple derivation of the 
expression in the main text for the 2-point function of the 
convergence, n. We start by assuming that the lensing 
occurs at late enough times that the anisotropic stress of 
the radiation can be neglected, so that in Newtonian gauge 
we can write the metric 



ds 2 = a 2 {i-j) [-(1 + 2§)drf + (1 - 2<S>)dx 2 } 



(Al) 



to first order in the gravitational potential $ ~ 1CP 5 . Here 
drj = a(t)dt is the conformal time and we have written the 
3-metric schematically as dx 2 — d\ 2 + r 2 (x)d£l. For scales 
smaller than the curvature scale we can approximate this 
as flat, r(\) = X, however on cosmological scales we need 
to use r = \K\~ X / 2 sinh li^ 1 / 2 ;^ for an open universe and 
r = \K\~ X I 2 sin \K \ l l 2 x for a closed universe. The con- 
formal factor, 0(77), accounts for the cosmological redshift 
of photon energy. When following null geodesies we may 
scale it out, i.e. set a = 1. The Lagrangian describing 
geodesic motion is L = ^g tlu x^x v where overdot represents 
differentiation with respect to an affine parameter A along 
the path. Recalling that the momentum p± = dL/dx± the 
Euler-Lagrange equations become 



dp± 
dX 



dL 

dx± 



= -2 



9$ dx 



dx 



-P\\ 



II 
d\ 



C($ 2 



(A2) 



Thus the deflection angle, a = Ap±/p^, receives differ- 
ential contribution da = — 2Vj_$c?a;||. Simple geometry 
dictates that such a deflection at a "lens" position results 
in a change of angle at the observer of 59 = [r LS /rs)5a 
where r^s = r(xs — Xl) is the (radial) distance from the 
lens to the source and r$ = r(xs) is the (radial) distance 
from the observer to the source. Translating this into a 
change in position at xl of 5x± = D^SO we see that two 
rays initially seperated by D$A.9 have a final seperation 



where 



Axt = (5ij - tjHj) DsA9j 



=2 j XS d x r -±^ didj® 

Jo r s 



(A3) 



(A4) 



The 2x2 matrix (Sij — ipij ) can be expanded in Pauli 
matrices with coefficients 



(s^ - i>ij) = (1 - k)i - 71 0-3 - 72^1 , 



(A5) 



so e.g. k — iTr(/i/>) = jipjj which leads to Eq. (||). Since 
all of the coefficients are derived from one function, spec- 
ifying any one of them is sufficient. We shall focus here 
on the convergence, k. Replacing with V 2 in the in- 
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so we can use 



tegral results in errors of £>(<!>) 
V 2 $ = 4irGpa 2 8 to obtain Eq. (g). 

To calculate the 2-point function of k we expand S(x) in 
Fourier modes and use the Rayleigh expansion of a plane 
wave to find 



C f = 4tt 



3 -n 

2 r - 



-1 2 



HA 



s(xi) sOca) 



«i 



«2 



y A ma SS ( fc ) J dXl 

ji{kxi)jt{kX2) 



dX2 



(A6) 



where g(x) is the distance kernel in Eq. ( JA4| ) . For open 
universes replace ji with the hyperspherical bessel func- 
tion. On small scales we may use the resolution of the 
identity 



k 2 dk je(kxi)je(kx2) 



^[Hx)]- 2 S(xi 



X2) 



to obtain the power per logarithmic interval in I, 



1(21 + l)C e _ 9?r 
4^ ~ ~i 



dx r 



9(X) 



i 2 



A 2 



(A7) 

Mr) ■ 
(A8) 
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In a spatially flat universe (r = x)i this reduces to Eq. (||). 
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